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We have developed a kinetic theory of hard needles undergoing binary collisions with loss of energy 
due to normal and tangential restitution. In addition, we have simulated many particle systems 
of granular hard needles. The theory, based on the assumption of a homogeneous cooling state, 
predicts that granular cooling of the needles proceeds in two stages: An exponential decay of the 
initial configuration to a state where translational and rotational energies take on a time independent 
ratio (not necessarily unity), followed by an algebraic decay of the total kinetic energy ~ t~ 2 . The 
simulations support the theory very well for low and moderate densities. For higher densities, we 
have observed the onset of the formation of clusters and shear bands. 



I. INTRODUCTION 

Gas kinetics of inelastically colliding particles has re- 
ceived a lot of interest in recent years, mainly in the con- 
text of granular matter |Q. Boltzmann's equation has 
been generalized to inelastic collisions which are char- 
acterized by normal restitution and possibly tangential 
friction or restitution Extensive simulations have 
been performed, in particular event driven (ED) algo- 
rithms are very effective for the kinetic gas regime. A 
variety of interesting phenomena have been observed, for 
example instabilities of the homogeneous state towards 
shearing or clustering . Most studies so far have con- 
centrated on spherically symmetric objects, whereas real 
grains arc in general nonsphcrical and often randomly 
shaped. The question arises whether the observed phe- 
nomena are generic for granular matter or model specific. 

In this paper we discuss the cooling properties of hard 
needles in terms of a time evolution operator, which ac- 
counts for the exchange of translational and rotational 
energy as well as for normal and tangential restitution. 
In addition we have performed ED simulations of large 
systems with up to 20000 needles. 

For low and moderate densities the system does not 
show any clustering instabilities but remains homoge- 
neous on the longest time scales, when the energy has 
decayed to 10~ n of its initial value. This allows us to 
formulate an approximate kinetic theory, based on the 
assumption of a homogeneous state. Cooling is found to 
proceed in two stages: 1) A fast exponential decay to a 
state which is characterized by a time independent ratio 
c of translational to rotational energy. 2) A slow alge- 
braic decay like t~ 2 of the total kinetic energy, while the 
ratio c exhibits small fluctuations around its time inde- 
pendent mean value. The latter is determined by the 
coefficient of restitution and by the distribution of mass 
along the needles, including equipartition for one par- 
ticular mass distribution. Simulations and approximate 
analytical theory are found to agree within a few percent. 

For high densities we observe large scale structures in 
the velocity field, similar to those seen in dense systems of 
smooth spheres &-@]. These structures give rise to long 



range correlations in the velocity field and the build up 
of large density fluctuations. No alignment of the needles 
with the hydrodynamic flow field is observed. 

II. DYNAMICS OF COLLISIONS 

Consider two rods of equal length L, mass m and mo- 
ment of inertia /. The center of mass coordinates are 
denoted by r± and r 2 . The orientations are specified by 
unit vectors U\ and u 2 , which span a plane E\2 with 
normal 



\U\ X u 2 \ 

We decompose r 12 = T\ — r 2 into a component per- 
pendicular r^ 2 — (ri2 ■ u±)u±_ and parallel r| 2 =: 
(S12W1 — S21U2) to E12 (see Fig. Q). The rods are in 
contact [0J|] if r± 2 = and simultaneously | s 12 1 < L/2 
and I sal I < L/2. 




FIG. 1. Configuration of two needles projected in the plane 
E12 spanned by the unitvectors iti and 112 

We want to determine the postcollisional center-of- 
mass velocities (v[, v' 2 ) and angular velocities (w'i, u> 2 ) 
in terms of the precollisional velocities (vi, V2, Wj., and 
u> 2 ) or momenta (p 1 = mvi,p 2 — mvz). Conservation 
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of total linear momentum implies 
p'i = Pi + Ap, 

p' 2 = p 2 - Ap. (2) 

No torque acts at the points of contact, so that 

= u>i + ^f-ui x Ap, 
J 2 =u> 2 - ?fu 2 x Ap. (3) 

holds. To determine the collision process we consider the 
relative velocity of the contact points which is given by 



Pi ~ P2 

V r = h Si 2 1ti - S 2 lU 2 . 



(4) 



We introduce two parameters e and (3 to characterize 
normal and tangential restitution: 

V' r ■ ux = -eV r -UX £€ [0, 1] (5) 

V' r ■ Ul = -pV r ■ Ul (6) 

V' r ■ u 2 = -PV T ■ u 2 6 [-1,1] (7) 

The above equations characterize the collision process 
completely and determine the postcollisional momenta 
in terms of the precollisional ones. 

In this paper we specialize to the case of perfectly 
smooth needles, i.e. j3 — — 1, which considerably simpli- 
fies the collision rules. The change of linear momentum 
is then given by Ap = au± with 



-(l + e)(V r -uj.) 
2/m+(s? 2 + sl 1 )// 



(8) 



III. ANALYTICAL THEORY 

The derivation of a pseudo-Liouville-operator for a 
system of N colliding, hard rods proceeds similar to the 
case of hard spheres [ptPPlOj] ■ The rods are confined to a 
3-dimensional volume V and interact via a hard core po- 
tential. The velocity of orientation ii\ = u>i X u\ is con- 
fined to the plane perpendicular to U\ and is therefore de- 
scribed by two generalized canonical momenta, pg 1 = I9\ 
and p ( f >1 = I(f>i sin 2 0\ , using spherical coordinates for the 
orientation. The total kinetic energy then reads 



N ( 1 1 

= E ( ^ + yA 

i=i y 
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(9) 



The time development of a dynamical variable A — 
A({r i (t),M i (t),Pj(t),'Ui(i)}) for positive times is deter- 
mined by the pseudo-Liouville-operator £ + 

A({ri, u t ,p^ in], t) = exp(i£+t)A({ri, m, p i: u z }, 0). 

(10) 

The pseudo-Liouville-operator consists of two parts, 
C + = Cq+C + . The first describes free streaming and can 



be expressed by the Poisson bracket with the kinetic part 
of the Hamiltonian iCo = {7ikin, • • • }p.b.- The second, 
£ + , describes binary collisions of two hard rods 



iC' 
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0(L/2 - \ Sij \)e(L/2 - \ Sji \)S(\ri]\ - + )(b t] - 1). (11) 

Here Q(x) is the Heaviside step function. The interpre- 
tation of the pseudo-Liouville-operator of eq. ([Ll]) is 
intuitively clear. The factor is the component 

of the relative velocity of the contact points perpendicu- 
lar to both rods. It yields the flux of incoming particles. 
®( — lt\ r ij\) 1S nonzero only if the two particles are ap- 
proaching and 0(L/2 - |fly|)0(L/2 - \sji\)5{\ri-\ - 0+) 
specifies the conditions for a collision to take place. The 
operator bij replaces momenta and angular momenta of 
particles i and j before collision by the corresponding 
ones after collision. 

The time evolution of nonequilibrium expectation val- 
ues of an observable A({r i ,u i ,p i ,u i },t) — A(T;t) is de- 
fined by: 

(A) t = J dTp(T;0)A(T;t)= J dTp(T;t)A(T; 0) . (12) 

r denotes the whole phasespace and p(T; t) is the N- 
particle phase space distribution function, whose time 
evolution p(T;t) = exp (— i£+t) p(T; 0) is governed by 
the adjoint C\_ of the time evolution operator C + . Here 
we are interested in the average translational and rota- 
tional kinetic energy per particle E tT = m/ {2N)^ ji v'l 
and E Iot = I/(2N) J^. u>f, as well as the total kinetic 
energy E = E tT + E Iot . 

It is impossible to calculate expectation values as given 
in eq. ([l2]) exactly and we are forced to approximate 
the N-particle distribution function. We assume that the 
system stays spatially homogeneous and that both lin- 
ear and angular momenta are normally distributed. In a 
system which is prepared in a thermal equilibrium state 
the initial decay rates can be computed exactly and yield 
different values for averaged translational and rotational 
energy. This suggests defining two different temperatures 
for the translational and rotational degrees of freedom, 
corresponding to the following ansatz for the N particle 
distribution function fl(i| 



PHCs(r;i) ~ exp 



E 



Ttr(t) T mt {t) 



(13) 



PHCs(r;^) depends on time via the average translational 
T tr (i) = 2/3(£ tr ) and rotational energy T mt (t) = (E mt ). 
We are interested in the cooling properties of a gas of 
hard needles and compute the expectation values T tT = 
2/3(i£ + E tI ) and T rot = (i£+E Iot ). Using the approxi- 
mate many particle distribution of eq. (|l3|), we find two 
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coupled differential equations 



V. HYDRODYNAMIC QUANTITIES 



2T tr 



I d 2 r^ ' T ' 



(l+T^kr 2 ) 1 ' 2 



1 + kr 2 



1 + e 



(1 + ^kr 2 f/ 2 
(1 +'fcr 2 ) 2 



, (14) 



4Trot - = - I d 2 r T - 



^kr 2 {l + ^kr 2 ) 1 ' 2 



3 7 T t 3 / 2 (l + e) 7n 1 + fcr 2 

l + e /• j5 fcr 2 (l + ^kr 2 ) 3 ' 2 



+ 



2 7n 



(l + £r 2 ) 2 



(15) 



with7= {2NL 2 yf¥)/{W^/m) and ft = (mL 2 )/(2I). The 
two dimensional integration extends over a square of unit 
length, centered at the origin. 



IV. SIMULATIONS 

Simulations are performed using an event driven algo- 
rithm where the particles follow an undisturbed trans- 
lational and rotational motion until a collision occurs. 
The velocities after the collision are computed according 
to the collision rules eqs. (^)-(^). For this algorithm col- 
lision times for each pair of needles have to be determined 
numerically, where we follow the algorithm proposed by 
Frenkel and Maguire 0. More efficiency is achieved by 
using the stratagem of Lubachevsky |ll[] and a linked cell 
structure, which allows us to look for collision partners 
only in the neighborhood. The algorithm is reasonably 
fast as long as there are only few needles in each cell 
of the linked cell structure, so that the time consuming 
search for collisions is restricted to the needles in the own 
and the neighboring cells. On the other hand we have to 
choose the linear dimension of these cells to be larger 
than the length of a needle, so that for high densities 
there are many needles in each cell and the algorithm 
becomes slow. 

We mention here that the algorithm runs into numer- 
ical problems if the time between two collisions becomes 
too short to be resolved properly as it usually happens 
during an inelastic collapse. To circumvent this problem, 
we use the i c -model |I2|] : if the time between a collision 
and the preceding one for at least one particle is smaller 
than a critical value t c , e is set to 1. We believe that the 
influence of this procedure is small: there occurred only 
two instances in the simulations presented here. 

We performed simulations for various system sizes 
(N = O(10 4 )) and coefficients of restitution in the regime 
of small and moderate densities (yL 3 < 1). For high 
densities (yL 3 > 10) only a few simulations could be 
done. We show here a simulation of N — 20000 needles 
in a box of length 12L with e = 0.9. 



To investigate deviations from homogeneity, we di- 
vide the simulation box into cells whose linear dimen- 
sion is chosen to be large compared to the length of 
the needles but small compared to the box size. Given 
the limitations due to finite size we choose cells such 
that on average about 25 needles are in one cell. We 
then compute for each cell the number density p a — 
X^ecciia 1 =: (!)<*) the translational energy per par- 
ticle PaE^ — (-Y v i)a and the hydrodynamic tempera- 
ture T* r = - mU 2 a /2. 

Momentum organization of the particles shows up 
in the hydrodynamic flow field p a U a — (vi) a . A 
good indicator |l3fl for the build up of long range ve- 
locity correlations is the ratio of the total kinetic en- 
ergy of the flow to the total internal energy S := 

(£ Q f^a)/(£ Q ^ r )- 

To investigate orientational ordering we compute in 

addition the quadrupolar moment of the needles Qf := 

(3(uj • e) 2 — 1). The unit vector e is chosen either along 

the direction of the particle's velocity or fixed in space. 



VI. SMALL AND MODERATE DENSITIES 

For densities such that yL 3 < 1, the system remains 
homogeneous, orientationally disorderd and without long 
range velocity correlations up to the longest observed 
time scales, i.e. when the energy has decayed to 10~ n 
of its initial value. To check for spatial clustering, his- 
tograms of fluctuations of the local density, velocity and 
translational energy were compared to those of an elastic 
system but no significant difference could be observed. 
Fluctuations, e.g. of the local density do not increase 
with time but remain stationary, so that we can compare 
our approximate theory with the simulations. We show 
here a simuation which has been performed for 10000 nee- 
dles in a box of length 24 L. This corresponds to a density 
of yL 3 w 0.72 or an average center of mass separation 

In this range of densities, cooling of a gas of hard nee- 
dles proceeds in two stages. First, there is an exponen- 
tially fast decay towards a state which is characterized 
by a constant ratio of translational and rotational energy. 
Second, there is an algebraically slow decay of both, the 
translational and rotational energy, such that their ratio 
remains constant. Both of these regimes are correctly 
predicted by our approximate theory, based on the as- 
sumption of spatial homogeneity. 

In Fig. |2| we plot the numerical solution of eqs. ( 14 Jl5|) 
for e = 0.8 and k = 6 as a function of dimcnsionlcss time 



r = 7 ^T tr (0)i. The total kinetic energy E = §T tr + T rot 
(in units of T tr (r = 0)) and the ratio T tr /T r 

ot are com- 
pared to simulations. Analytical theory and simulation 
are found to agree within a few percent over eight orders 
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of magnitude in time, 
initial condition). 



(2rot(0) = has been chosen as 




FIG. 2. Double logarithmic plot of total kinetic energy 
E = 3/2T tr + T rot in units of T tt (r = 0) and of T tI /T rot 
versus dimensionless time r. The simulation data are from a 
system of 10000 needles, box of length 24 L and e = 0.8. 

The decay of T tr /T rot to a constant value happens on 
a timescale of order one. In this range of times the total 
kinetic energy E remains approximately constant (on a 
logarithmic scale) and decays like t~ 2 only after trans- 
lational and rotational energy have reached a constant 
ratio. 

Eqs. (|l4|,|l5]) allow for a solution with a constant ratio 
of Tt r /T r ot = c and both T tr and T Iot decaying like i~~ 2 . 
To determine the constant c, we plug the ansatz cT rot = 
T tr into eqs. (pp5|) and use cT rot - f tT = 0. This yields 
an implicit equation for c, whose solution is plotted in 
Fig. |3| as a function of k and e. 
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FIG. 3. Asymptotic ratio Ttr/Trot as a function of e and k. 
Setting c = 1 in this implicit equation yields an equa- 



tion for k which reads 
(1-e 2 ) 



2^ / d 2. 



_ 1 - ^k*- 2 
\/l + k*r 2 



= 0, 



(16) 



i.e. equipartition holds for all values of e < 1 if k = 
(mL 2 )/(2I) is set to the particular value k* = 4.3607, 
given as the solution of eq. ([l6|). For e = 1, equipartition 
always holds, independent of k. 

For k < k* we find T tr < T rot and for fc > k* T tr > T rot . 
Hence the distribution of mass along the rods determines 
the asymptotic ratio of rotational and translational en- 
ergy, including equipartition as a special case. 

The asymptotic solution discussed above is approached 
for arbitrary initial conditions for long times. If a totally 
elastic system is prepared in an initial condition with 
T tr 7^ T rot , we expect that the equilibrium state (equipar- 
tition) is reached exponentially fast with a relaxation rate 
given by v ~ jy/E(0). As long as energy dissipation due 
to inelastic collisions is small, we expect similar behavior, 
as indeed the numerical simulations show. 



VII. DENSE SYSTEMS 

Simulations of inelastic hard spheres show well devel- 
oped density clusters and vortex patterns ^-gQ, if al- 
lowed to evolve freely for sufficiently long times. The 
dominant mechanism for the formation of vortex struc- 
tures has been traced back to noise reduction [|l5| : After 
many collisions the particles move more and more par- 
allel. It is not clear a priori, whether such a mechanism 
should also apply to rotating needles. In the simulations 
we clearly observe the formation of large scale structures 
in the velocity field. In Fig. ^ we show the hydrdynamic 
flow field after 600 collisions per particle for a system of 
density ^-L 3 w 11.6, corresponding to an average center 
of mass separation l/L = 0.44. 
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FIG. 4. Flow field of the system (20000 needles in a box of 
volume (12L) 3 and e = 0.9) after 600 collisions per particle 
(the length of the velocity vectors are in arbitrary units). 

We observe two shear bands (note the periodic bound- 
ary conditions), which move in opposite directions. 
Within a band the local flow field is to a large degree 
aligned. The dominant part of the velocity of each par- 
ticle Vi is given by the flow U so that a large fraction of 
the kinetic energy is in the flow and the ratio S should be 
high. In Fig. |5| a) we show S as a function of the number 
of collisions per particle. S increases from a value of 0.05 
to a value of about 2.5, i.e. by a factor of 50. 
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collisions/particle x 

FIG. 5. a) S as a function of collisions per particle, b) 
Double logarithmic plot of kinetic energy per degree of free- 
dom T tr and T Iot in units of T tr (r = 0) versus dimensionless 
time r. 

To visualize spatial inhomogeneities we plot in Fig. || 
regions with local deviation of the density Ap a > 0.5 at 
the beginning of the simulation and after 600 collisions 
per particle. Obviously clustering occurs. To quantify 
this obseravtion we have computed the second moment 
of the density fluctuations. It is found to increase by 



a factor of 6 over its initial value after 400 collisions. 
After about 500 collisions per particle it decreases again, 
indicating that clusters form and dissolve again. 
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FIG. 6. a) Density at the beginning of the simulations. 
Only very few and small regions have a more than 50 % higher 
density than average, b) Large regions with higher density 
have built up after 600 collisions per particle. 



For spheres one observes 14 1 that most of the mass 
is concentrated in the two counterflowing streams. To 
check for correlations between flow field and mass density, 
we plot in Fig. ^ the components of the flow field U y , U z 
and the density fluctuation as a function of x, for fixed 
y = 6L and averaged over 10 values of z = 1.2L, 12L. 
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This and similar plots give no hint of a strong correlation 
between flow field and density fluctuations. 




FIG. 7. y— and z— component of the flow field and fluc- 
tuation of the density as a function of x, for fixed y — 6L, 
averaged over 10 z- values. 

We also observe a weak tendency towards organisation 
of rotational velocities. The ratio of the kinetic energy of 
the local rotational flow to the local rotational temper- 
atur is found to increase by about 50 % (as compared to 
an increase by a factor of 50 for the translational veloc- 
ity). Consequently, the deviation from Haff's t~ 2 cooling 
law is much stronger for the translational degrees of free- 
dom T t r than for the rotational degrees of freedom T rot 
(see Fig. |b)). 

To investigate alignment of the needles with the large 
scale velocity flow field, we compute the quadrupolar mo- 
mentum Qf with respect to the particle velocity, i.e. 
e = Vi/\vi\. A histogram over all needles is shown in 
Fig. H The configuration after 600 collisions per parti- 
cle is compared to the initial state which corresponds to 
randomized orientations. In addition we plot the theoret- 
ical prediction for the histogram (straight line) which has 
been calculated on the assumption that rods are oriented 
randomly and independent of their velocity. No indica- 
tion for alignment of the needles can be seen. Neither do 
we observe a tendency for global ordering. 
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FIG. 8. Histogram of Q Vi . The distribution after 600 col- 
lisions per particle coincides with the initial distribution and 
with the prediction for randomly distributed orientations. 



VIII. CONCLUSION AND OUTLOOK 

Our aim was to systematically investigate the cool- 
ing dynamics of a granular system of nonspherical ob- 
jects. We have chosen the simplest nonspherical grains, 
hard needles, for which we were able to formulate an 
approximate kinetic theory, generalising methods of ki- 
netic theory of elastic systems to inelastic collisions 
with normal and tangential restitution. In addition, sim- 
ulations of large systems for various densities have been 
performed. 

Analytical theory is so far restricted to the regime 
of moderate densities, where the interparticle spacing 
is comparable to or larger than the length of the nee- 
dles. Such systems remain homogeneous and show no 
long range correlations in the velocity field or orienta- 
tion. Cooling proceeds in two stages: 1) An exponentially 
fast initial decay towards a state with constant ratio of 
translational to rotational energy and 2) an algebraically 
slow decay, such that the above ratio remains constant 
in time. The ratio of translational to rotational energy 
is controlled by the coefficient of normal restitution and 
by the distribution of mass along the rods. 

Simulations in the dense regime, where the interpar- 
ticle spacing is smaller than the length of the needles, 
reveal large scale structures in the translational velocity 
field. The density does not remain homogeneous, but 
clusters form and dissolve again. Cooling proceeds in 
three stages. For short and intermediate time scales the 
relaxation is similar to the low density system, whereas 
on the longest time scales we observe a crossover from 
the algebraic decay to an even slower decay. This latter 
decay may be identical to the asymptotic cooling law of 
hard inelastic spheres, r~ d//2 in d— dimensions, where r 
is the average number of collisions suffered by a particle 
within time t, which has been derived in |13|. 
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We plan to generalise the hydrodynamic analysis to 
grains with rotational degrees of freedom and in partic- 
ular hard rods. Another possible extension of our work 
are rods of finite width and with spherical endcaps. 
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